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The interstellar medium (ISM) provides a unique laboratory for highly supersonic, driven hydro- 
dynamics turbulence. We present a theory of such turbulence, confirm it by numerical simulations, 
and use the results to explain observational properties of interstellar molecular clouds, the regions 
where stars are born. 



1. Introduction. Stars are formed as a result of grav- 
itation (Jeans) collapse of dense clumps in interstellar 
molecular clouds. The structure of such clouds in a 
large interval of scales (from about lOOpc to O.Olpc) lacks 
any characteristic length, and can be understood as aris- 
ing from supersonic hydrodynamic motions sustained on 
large scales by supernovae explosions A fluid 

motion with characteristic large-scale relative velocities 
of order 1-10 km/sec compresses rapidly radiating and 
therefore relatively cold molecular gas (T ^ IQK) up 
to very high densities (above lO^cm"'^). Instabilities of 
shock fronts create a hierarchy of gas clumps with broad 
mass, size, and velocity distribution controlled by the 
Mach number (M = v/c) and by the Alfvenic Mach 
number [Ma = v/va), where Va is the Alfven veloc- 
ity, Va = Bj^Amp, and c is the sound speed. Depend- 
ing on parameters of clouds, the Mach number can be 
about 30 on the largest scales, and the Alfvenic Mach 
number can exceed 1 as well [|o|j2|,|l|,|l|,||,||,|5|j2|l . 
A systematic study of scaling properties of supersonic 
turbulence and of the structure of molecular clouds was 
initiated by the work of Larson [^,^ , but despite the 
large number of observational and numerical investiga- 
tions that appeared for the last 20 years, the theoretical 
understanding of the turbulence has been rather poor. 

In the present paper we provide a new analytical model 
for a supersonic turbulent cascade and test the model 
against observations and numerical simulations. We find 
a very good agreement within the error bar uncertainty. 
The analytical approach is suggested by the following two 
numerical results. First of all, we establish that for large 
Mach numbers (M > 2) the velocity field in the inertial 
interval is mostly divergence-free and shear-dominated, 
with the intensity of its potential component being only 
about 10% of the intensity of its solenoidal part. The sec- 
ond finding is that the most intense dissipative structures 
of the turbulence are two-dimensional sheets or shocks 



as oppose to the incompressible case where most of en- 
ergy is dissipated in filaments, see also [^,^. Using 
these two ingredients in the framework of the so-called 
She-Leveque model of turbulence we calculate two-point 
correlators of velocity and density fields. This allows us 
to construct the multifractal distributions of velocity and 
density fields, that statistically describe the structure of 
a turbulent molecular cloud. In the present paper we 
present the analytic derivation of the model and sum- 
marize the most important numerical results. The de- 
tailed numerical and observational analysis will appear 
elsewhere P,p7t. 

In the following section we construct the multifractal 
distribution of the velocity field, and numerically check 
the first 10 velocity-difference structure functions [the 
definition is given below]. In section 3 we proceed with 
the density distribution, derive a general formula for two- 
point density correlators and compare the result with the 
numerical simulations. Conclusions, applications, and fu- 
ture research are outlined in section 4. 

2. Multifractal model of supersonic turbulence. In 1994 
She and Leveque suggested a model that turned out to 
be very successful in explaining the experimental find- 
ings for incompressible turbulence p^ . The model rep- 
resents a turbulent cascade as an infinitely-divisible log- 
Poisson process [p| j3^ , p^ , and has three input parame- 
ters. The first two of them are the so-called naive, i.e., 
non-intermittent, scaling exponents of velocity fluctua- 
tion and of the eddy turnover time. The non-intermittent 
velocity of the eddy of size I scales as vi ~ Z®, and the 
"eddy turnover" time scales as i/ ^ l^. The third pa- 
rameter is the dimension of the most intense dissipa- 
tive structure, D. In the incompressible case, the first 
two parameters take their Kolmogorov values, A = 2/3, 
Q = 1/3, and the third parameter is D = 1 since the 
dissipation mostly occurs in elongated filaments. 

The model predicts the so-called structure functions of 
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the velocity field that are defined as follows 

Sn{l) = {Hx + I) - V{X)\) ^ l^^^l (1) 

The velocity components in this formula can be either 
parallel or perpendicular to vector I. In the former case 
the structure function is called longitudinal, in the lat- 
ter case transversal. There is experimental evidence that 
both scale in the same way [|||Jl^, therefore, we will 
not distinguish between them as far as the scaling is 
concerned. The She-Levequc formula gives the follow- 
ing expression for the scaling exponents of the structure 
functions, 



C(n) = 9(1 - A)n + (3 - D){1 - 



(2) 



where S = 1 — A/(3 — Z?) Q]. This expression has been 
experimentally checked to work for structure functions 
up to the 10th order, within an accuracy of a few per- 
cent |,||. 

To address the supersonic case, we note that in the 
inertial interval the turbulence is still mostly divergence- 
free (Fig. (|l|)). The effect of vorticity generation in the 
random 3D flow is analogous to the effect of magnetic dy- 
namo, since magnetic field and vorticity obey the same 
dynamic equation. We therefore leave the Kolmogorov 
parameters Q and A unchanged. 
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FIG. 1. Results from numerical simulations of ran- 
domly driven MHD equations with resolution 128"^. The 
ratios of the potential, Ec, to the solenoidal,i5i, part of 
the velocity field are presented for the Mach numbers up 
to A/ — 10 at the largest scale. The initial magnetic Mach 
number has been chosen in the range Ma ~ 0.6, . . . , 10, 
and the numerical integration has been conducted for 
several turn-over times of the largest eddies, which was 
enough to reach the steady state. The isothermal equa- 
tion of state was used and the large-scale driving force was 
Gaussian, solenoidal, and correlated at a turn-over time 
of the largest eddy. The detailed description of the nu- 
merical set up can be found in [4]. The scaling relations 
reported in the present paper were obtained for M ~ 10, 
and Ma ~ 3. 



However, the most dissipative structures of a super- 
sonic flow are different from the incompressible case. 
Instead of filaments they look like shocks or two- 
dimensional dissipative sheets, therefore, D = 2 [^,1). 
With such input parameters, formula ^ is recast as 



Cin) =n/9 + l-(l/3) 



n/3 



(3) 



Quite remarkably, this formula works with good accuracy 
for the numerical simulations, fitting the first ten struc- 
ture functions with an error of about 5%, see Fig.(0). 



FIG. 2. Slopes of transversal structure functions com- 
puted for n = 1, 2, . . . , 10 (correspondingly, from bottom 
to top), in the 250^ run with M = 10 and Ma = 3. The 
plot presents the ratios of the differential slopes of the 
structure functions to the differential slope of S3{1). These 
ratios, ({n)/({3), exhibit excellent scalings, in agreement 
with the Extended Self-Similarity hvpothesis [2]. They 
are well described by our formula (W). Note the strong 
difference of the scalings of the structure functions from 
the scalings given by the Kolmogorov model, C{n) = n/3, 
and by the Burgers model, ("(n) = 1, see [15]. 

The observational data consistently indicate steeper 
than Kolmogorov velocity spectra, see e.g. [ESLH. Our 
model (|^) predicts jw^P fc~^~^(^) — fc^~^% which 
agrees with observations within error bars. Interestingly 
enough, the average velocity spectrum, originally inferred 
from observations by Larson, was k~^-'^^ on scales of 
order Ipc < I < lOOOpc, and k~^ '^^ towards smaller 
scales, O.lpc < I < lOOpc [|l|,|o). 

On the analytical side, the model (|^) implies the so- 
called multi-fractal distribution of the turbulent fluctua- 
tions. Here we discuss the statistics of the velocity field, 
and in the next section, apply the results to the den- 
sity field. To visualize the model, assume that the whole 
space (a molecular cloud or a simulation domain) con- 
tains turbulent structures of different (in general, fractal) 
dimensions. In the vicinity of a particular structure, the 
velocity difference scales with some particular exponent 
that we denote by h, i.e., vi ^ l^. The dimension of this 
fractal structure will be denoted D{h). If we divide the 
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space into small boxes of size /, then the number of boxes 
covering the fractal structure of dimension D is propor- 
tional to l"^ , while the total number of boxes is propor- 
tional to l^'^. The probability to find ourselves inside a 
box covering the fractal with dimension D{h) is there- 
fore Ph{l) ^ l3-D{h) ^ average the n's moment of the 
velocity difference we just need to sum F''^ph{l), which is 
the contribution of one particular fractal structure, over 
all the fractal structures. We get 

5„(0~^r''+'^^('')'-/^("^ (4) 

h 

Knowing D{h) is equivalent to knowing CI"-); these 
two functions are related by the Legendre transform as 
one immediately gets by assuming that I is small, and 
by evaluating the sum in (||) by the steepest descent 
method For example, knowing from (0), one 

can restore D{h); we, however, will not need D(h) for 
our present purposes. 

3. Multifractal distribution of the density field. First 
we would like to derive an important relation of super- 
sonic turbulence. We start with the Navier-Stokes and 
continuity equations: 

5tu+(u- V)u= (ry/p)Au-c2V/9/p + f, (5) 
dtp + W ■ (pu) ^ 0. (6) 

Let us introduce the density correlator R{t, x) — 
{p{xi,t)p{x2,t)), and the density-weighted second-order 
velocity structure function, G^'^{t,x) = {p{l) p{2)[u{l) — 
u{2)Y[u{l) — w(2)]'^), where x = xi — X2. Averaging is 
performed over the random force. Differentiating R{t, x) 
with respect to time, and using (||) and (^), one gets: 

d^R + V.VfcG'" - c^AR + V,([/(l) - /(2)]V(l)p(2)) 

= 277AV,(u(l)V(2)), (7) 

where the spatial derivatives are taken with respect to x. 
In the inertial interval, the forcing and the viscous terms 
are small. In the supersonic regime, one can also neglect 
the term. Assuming now that the turbulence is in 
a steady state, we are left with a simple equation that 
must hold in the inertial interval, ViVfeG**'' = 0. Due to 
spatial isotropy we get: 

G'''{x)=AS''' + 0{xyL^) + ..., (8) 

where L is the external force correlation length, and A 
is some constant. In the inertial interval, x <ti L, one 
gets G — const. To obtain the density distribution let 
us make a natural assumption that both the density and 
the velocity fields have fixed scalings in the vicinity of 
the same turbulent structures. In other words, close to 
a fractal structure where the velocity field has scaling h, 
the density field has some other, but also constant along 
the same structure, scaling i.e., p{l) P^''-'. The 

condition (||) now reads 



G^Y. - const. (9) 

h 

The other restriction comes from the mass conservation 
law, 

(p) ^^rW+3-^(M = const. (10) 
h 

Strictly speaking, our constraint conditions and ( p^ ) 
are to a certain extent phenomenological. However, we 
found them to be consistent with our simulations. More- 
over, the theory based on them predicts density corre- 
lators rather successfully, which we are going to demon- 
strate now. Let us assume that the function a{h) is ana- 
lytic and can be expanded as a{h) — a+bh+gh^+. ... As 
a minimal model, consider the case, when a{h) is a linear 
function, a{h) — a + bh. It turns out that such a linear 
ansatz is consistent with both restriction conditions (j^) 
and jiol). As follows from (^, the mass conservation con- 
ditionWo) is satisfied with a — ~C{p) and b = p, where p 
is arbitrary. The equation for p is then derived from the 
dynamic constraint (||), which gives C(2p -I- 2) = ^({p). 
The solution of this equation will be denoted as pq. If 
we use our formula (^, po can be found exactly; we thus 
obtain b ^ pa = 2.28 and a = -C(Po) = -0.82. The 
multifractal distribution for the density field, Dp{a), is 
thus related to the multifractal distribution of the veloc- 
ity field, Dp{a) = D [{a + C{po))/po]. Fractal and mul- 
tifractal distributions of density fields have indeed been 
inferred from observations, see, e.g., |p^,[7|, and from nu- 
merical simulations 

By analogy with the velocity field, the quantities of 
practical interest are density correlators. They can be 
calculated with the aid of a formula analogous to (||), 

{[pix + l)pix)]"') - J2 ^ ;«("0. (11) 

h 

Upon substituting the linear expression for a{h) and us- 
ing formula (Q), one immediately gets £,{m) = C(2mpo) ~ 
2rn^(po)- This formula allows one to obtain the density 
scaling if the velocity structure functions are known ei- 
ther from theory or from experiment. Since C(^) is a 
concave function, ^(m) is negative for m > 1/2. For 
TO = 1,2,3 the formula gives ^(1) ~ -0.3, ^(2) =i -1.3, 
and ^(3) ~ —2.4, values close to what is obtained in the 
numerics, see Fig. (||). We give here only the first three 
exponents since starting from to = 3, the density ex- 
ponents depend on velocity structure functions of order 
higher than 13, which cannot be reliably produced with 
our numerical resolution. 

4. Conclusions. We have suggested a self-consistent 
model that provides an explanation for numerical and 
observational scaling laws of supersonic ISM turbulence, 
the so-called Larson's laws. We would like to conclude 
with the following remarks: 
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1. The She-Leveque approach was also appUed to in- 
compressible MHD turbulence in [p^ , pl]j2^ ], where 
different scalings were suggested. Most success- 
ful was the approach of Miiller and Biskamp , 
where the energy cascade was assumed to be 
Kolmogorov-like, but the dissipation occurred in 
micro current sheets. In this case, the same for- 
mula (|^) gave a good agreement with numeri- 
cal simulations for structure functions up to or- 
der 8. Our results together with those by Miiller 
and Biskamp support the ideas put forward in ||] 
and ll^, that completely different turbulent sys- 
tems can belong to the same class of universality, 
i.e. have the same velocity scaling exponents. 

2. Turbulence with small pressure is usually re- 
ferred to as Burgers turbulence, the theory of 
which has been rapidly developing in recent 
years [|2|j37||Jl^,|l^,||,^|j3^ . However, the Burg- 
ers velocity field is usually assumed to be poten- 
tial, which is true in one and two dimensions, 
but inconsistent with the 3D case due to strong 
vorticity generation. However, our general rela- 
tion, p(l)p(2) ~ \v{l) - w(2)|2po/Z2C(Po)^ which is 
valid inside any correlation function, can be useful 
for the closure problems of Burgers turbulence. 

3. Our model is consistent with available observa- 
tional results, although the error bars of observed 
velocity structure functions are too large for a pre- 
cise comparison. Moreover, only projected quanti- 
ties (i.e., integrated along the line of sight) are ob- 
servationally available, and therefore the 3D results 
should be reformulated for these projected fields — 
this is a subject of future work |27[. 
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FIG. 3. Density correlators for m = 1,2,3. The nu- 
merically obtained slopes are ^(1) ~ —0.3, ^(2) ~ —1.3, 
and ^(3) ~ —2.1, close to the theoretical prediction (pj|). 
The numerical simulations are the same as in Fig. {nh. 
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